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Abstract 

We describe a computational investigation of tunneling at finite 
energy in a weakly coupled quantum mechanical system with two de- 
grees of freedom. We compare a full quantum mechanical analysis to 
the results obtained by making use of a semiclassical technique de- 
veloped in the context of instanton-like transitions in quantum field 
theory. This latter technique is based on an analytic continuation of 
the degrees of freedom into a complex phase space, and the simultane- 
ous analytic continuation of the equations of motion into the complex 
time plane. 

1 Introduction and Motivation 

The existence of a small parameter ("coupling constant") in quantum me- 
chanical systems leads to a (typically asymptotic) expansion of observables 
in powers of this parameter. No other approximation technique has proved 
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as powerful in obtaining physical predictions in such diverse fields as atomic 
physics, chemistry, quantum field theory, etc. Despite these successes, many 
phenomena in such systems are not amenable to perturbation theory: for 
example barrier penetration in quantum mechanics does not occur at any 
order in an (asymptotic) expansion in powers of h. 

Techniques for dealing with non-perturbative phenomena in theories with 
a small parameter are far less general. Perhaps the best known example is 
the WKB approximation, familiar from one-dimensional wave mechanics. 
A similar technique, the instanton method, is often used to discuss certain 
non-perturbative phenomena in quantum field theory. But in more compli- 
cated cases these methods fail: examples are tunneling at non-zero energy 
in quantum mechanical systems with more than one degree of freedom, and 
tunneling processes in quantum field-theoretic models with exclusive initial 
or final states. 

Recently the authors of Ref. |], [| have suggested a method for dealing 
with high-energy processes that proceed through tunneling in weakly coupled 
quantum field theory. Their technique begins with a path integral represen- 
tation of the matrix element in question, followed by a double analytic con- 
tinuation: the fields in the path integral are continued to the complex plane, 
and in addition the time evolution is continued along a complex contour. In 
spite of these complications, this technique is essentially semiclassical. The 
resulting complexified classical system typically remains intractable to an- 
alytic methods. Consequently computational techniques must be employed 
to obtain quantitative results; the feasibility of the corresponding calcula- 
tions in field theory has been demonstrated in Ref. [[| [|. We should also 
stress that the validity of the formalism of Ref. [JTJ, [| has not been proven, 
though its plausibility has been supported by comparison with perturbative 
calculations about the instanton |(J . 

More specifically, the process under discussion is a non-perturbative in- 
stanton mediated transition induced by the collision of two highly energetic 
particles. Perturbative calculations (see for instance Ref. f?|, ||, |j and refer- 
ences therein) about the instanton suggest that the total cross section has 
the following functional form 

a 2 ^ any oc e ? (1) 



where g is the small coupling constant of the theory and E is the center- 
of-mass energy. To compute the leading exponent the authors of Ref. [2|] 



suggested considering an inclusive process with a large number n of incoming 
particles. They argued that the total probability has a similar form 

CTn^any OC t s* (2) 

and that the exponent F(g 2 E, g 2 n) can be calculated semiclassically by con- 
sidering a complexified classical system. Furthermore, they conjectured that 
the two-particle exponent Fq in Eq. (|1|) is an appropriate limit of the multi- 
particle one: 

F (g 2 E) = lim g 2 n ^ F(g 2 E,g 2 n) (3) 

Equations ([!]), (§) on the one hand, and equation (§) on the other, have 
different status. While the validity of Eq. (0) has been demonstrated by 
path integral methods (cfr. Section 4), neither the general functional form 
([]]) nor the limiting procedure (Q) have been proven so far. 

Since the formalism of Ref. |], § has not been rigorously derived from 
first principles, and the direct evaluation of the resulting path integral, by 
computer simulation or other numerical procedures is beyond current reach, 
we have chosen to test the technique by reducing the number of degrees of 
freedom. In quantum field theoretic models, a general field configuration 
may be expanded in a complete (infinite) basis of normal modes. In the 
asymptotic time domains t — > ±oo these modes are non-interacting, and the 
evolution is characterized by a definite particle number. In a semiclassical 
description of the tunneling process the field evolves through a non-linear 
regime, and the crucial question is how the particle numbers in the incoming 
and outgoing asympotic states are related by the non-linear evolution. A 
minimal model capable of mimicking this dynamics will have some internal 
degree of freedom, whose excitations at asymptotic times will correspond to 
the particle number of the field theoretical system, and a non-linear inter- 
action with a barrier, that can be penetrated by tunneling. This can be 
realized with a system of two particles moving in one dimension. Let the 
coordinates of these particles be x\ and 22, and the dynamics be described 
by the Lagrangian: 

L = ^xl + - l -u?{x x - x 2 f - V{ Xl ) (4) 

where V is an arbitrary positive semi-definite potential which vanishes asymp- 
toticallyQ Since the theory is to be weakly coupled, we assume a potential 

1 We could of course allow V to depend on xi as well, provided it does not depend only 
on the combination x\ — Xi. 



of the form 



V(x) = \u(gx) 



(5) 



with gCl. For simplicity we will use a gaussian for the potential 



U(x) = e 2 



(6) 



although the treatment of other potentials is similar. The properties of the 
system described by the above Lagrangian are made clearer by replacing the 
variables xi,x% with the center of mass coordinate X = (x± + x^)/2 and the 
relative coordinate y = (x\ — x<i)j1. With this substitution the Lagrangian 
takes the form 



and we see that asymptotically it describes the free motion of the center of 
mass and a decoupled harmonic oscillator. Within the range of the potential, 
though, the two degrees of freedom are coupled, giving rise to a transfer of 
energy between them. 

In the classical case, the coupling g is an irrelevant parameter: we may 
rescale the degrees of freedom so that g appears as a universal multiplicative 
factor. Defining new coordinates X = gX, y = gy the Lagrangian becomes 



The value of g is crucial, however, for the quantum system: the path integral 
formulation of quantum mechanics together with Eq. @ show that g 2 plays 
a role similar to h in determining the magnitude of the quantum fluctuations; 
the classical limit corresponds to g 2 — > 0. This is in close analogy with the 
field theoretical systems mentioned above. In the following we will use units 
with h = 1 and will characterize the semiclassical treatment as an expansion 
for small g. 

The repulsive potential implies a barrier that must be either overcome 
or penetrated through tunneling for a transition from an initial state where 
the center of mass coordinate is approaching the barrier from, e.g., large 
negative X, to a final state where it is moving away from it towards large 
positive X. The corresponding transmission probability T will depend on 
three quantities: the total initial energy E\ the initial energy of the oscillator 
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E OS c (or, equivalently, on its initial quantum number n related to its energy 
by E osc = (n + l/2)u); and the value of g. In analogy with Eq. (JJ), the 
transmission probability for an oscillator initially in its ground state as g — > 
has the asymptotic form 

%(E) = C (g 2 E)e-^ Fo{92E) (9) 

for some prefactor Cq. Likewise, for a transition from an initial state in the 
n-th excited level one expects, for g — > 0, ng 2 fixed, 

T n (E)=C(g 2 E,g 2 n)e-^ F{92E ' 92n) (10) 

in analogy to Eq. (^). The advantage of the model we are considering is 
that the process in question admits a full quantum mechanical treatment as 
well as the semiclassical analysis. We will present the results of a numerical 
solution of the full Schrodinger equation and show that the transition from 
the oscillator ground state can indeed be fitted very well with the expression 
of Eq. (H). Independently of the semiclassical analysis, this result represents 
a direct verification of the functional form of Eq. (|I]), (^). This will be our 
first conclusion. 

We will then use the technique of Ref. |l|, H and evaluate the function 
E(g 2 E, g 2 n) entering Eq. (|1(]) by solving numerically the complexified classi- 
cal equations on the appropriate contour in the complex time plane. We will 
thus be able to check the validity of Eq. (|3|), with the l.h.s., F (g 2 E), obtained 
through the full quantum mechanical treatment and the r.h.s., F(g 2 E, g 2 n), 
calculated in a semiclassical way. We will show that Eq. (^) indeed holds, 
and so, in the context of our model at least, we will be able to confirm the 
conjecture of Ref. JL], U by a direct numerical computation. This will be the 
second main conclusion of this paper. 

2 The Classical System 

Let us first consider a classical evolution whereby the two particles are ini- 
tially located on the negative x-axis well outside the range of the potential 
and their center of mass is moving with positive velocity (i.e. toward the bar- 
rier). The motion of the system is specified completely by four initial value 
data. Time translation invariance of the system allows us to choose one of 
these to be the initial time. It is convenient to take the remaining three 



to be the rescaled total energy of the system, e = g 2 E, the rescaled initial 
oscillator excitation number, v = g 2 n, (in the classical theory n is defined 
as E osc /u and need not be integral) and an initial oscillator phase, <fi. The 
question at this stage is whether the system can cross to the other side of the 
barrier, i.e. whether the transition is classically allowed. In particular, in the 
projection to the v-e plane there will be a classically allowed region where, for 
some value (s) of 0, the system will evolve to the other side of the potential 
barrier. The rest of the plane will consitute the classically forbidden region 
where, no matter what the initial phase, the system will bounce back from 
the barrier. 

Clearly the entire domain e < 1 belongs to the classically forbidden region: 
there can be no classical transition with a total energy smaller than the 
barrier height (equal to 1 in rescaled units). However, a total energy larger 
than the barrier height is per se no guarantee that the system will cross to 
the other side of the barrier. The coupling between the center of mass and 
oscillator degrees of freedom due to the potential will cause a transfer of 
energy between the two, whose net effect can be repulsion from the barrier 
even when the total energy is larger than the barrier height. In general, for 
every initial value of v there will be some minimal rescaled energy eo(^) such 
that for e > eo transitions across the barrier are possible. The function eo(^) 
describes the boundary of the classically allowed region. 

The minimum of e^iv) is equal to 1 (i.e. to the barrier height). Indeed, 
there is an obvious, unstable, static solution of the equations of motion with 
both particles on top of the potential barrier (xi(t) = X2(t) = 0). This solu- 
tion, incidentally, corresponds to the static solution called the "sphaleron" in 
instanton mediated processes [f[I]] . If one perturbs this solution by giving an 
arbitrarily small, common positive velocity to both particles, they will move 
in the positive direction toward X = oo. (It is easy to prove that the parti- 
cles cannot go back over the barrier in this situation. If this were to happen, 
at some moment in time x\ would pass through zero. At that moment, by 
conservation of energy, the magnitude of the center-of-mass velocity could 
not be larger than the initial velocity of the particles. But a perturbative 
analysis of the initial motion shows that, however small its initial velocity 
may be, the center of mass will acquire some finite positive momentum, which 
will continue to increase so long as x\ > 0. This implies that the magnitude 
of the center of mass velocity cannot revert to its original arbitrarily small 
value.) Similarly, the time reversed evolution has the two particles proceed- 
ing towards X = — oo. The two evolutions, combined, describe therefore a 



classical process where the system goes over the barrier with an energy larger, 
but arbitrarily close to the barrier height. This evolution, obtained by an 
infinitesimal perturbation of the "sphaleron" , will produce a definite asymp- 
totic value uq of the rescaled initial excitation number, which will characterize 
the minimum of €q(v). 

2 



1.5 



u 1 



0.5 



"0 12 3 4 

v 

Figure 1: Boundary of the region of classically allowed transitions. 

Values of v smaller or larger than uq will give rise to values eo(u) larger 
than the barrier height. Of particular interest for us is e (0), i.e. the lowest 
energy for which one can have a classically allowed transition with the in- 
coming system in its ground state. In the limiting cases of infinite oscillator 
strength (uj — > oo) or of decoupled particles (u = 0) eo(0) takes values 1 
and 2, respectively. Indeed, the former case is equivalent to having a single 
particle (with twice the mass), which will evolve over the barrier as soon as 
its initial energy is larger than the barrier height. In the second case, the 
particles proceed independently, sharing the initial energy, and the center 
of mass will move to positive infinity whenever the particle that feels the 
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potential will be able to move over the barrier, which of course will happen 
when its share of the total energy is larger than 1. For finite, non-vanishing 
values of lu, the value of eo(0), or, more generally, e (z/), must be determined 
numerically The small number of degrees of freedom in our model allows us 
to explore the phase space systematically, varying <p, e and v. In this way we 
can determine the boundary of the classically allowed region with reasonable 
accuracy. Throughout this paper we will take u = 1/2. The corresponding 
boundary of the classically allowed region is illustrated in Fig. p]. (The dotted 
line in the figure represents the kinematic boundary ujv < e.) With u = 1/2, 
eo(0) equals approximately 1.8. For energies lower than this value, an incom- 
ing system in its ground state may transit the barrier only through tunneling. 
In the next section we will present a full quantum-mechanical calculation of 
the corresponding transmission probability Tq(E, g). In Sect. [| we will adapt 
the semiclassical formalism of Ref. [fll to derive an exponential expression 
for the leading factor in To, and calculate the exponent. 

It is worthwile to mention that the existence of a definite eo(0) where the 
boundary of classically allowed region meets the axis v = is one important 
property which our simplified model, most likely, does not share with its 
more complex field theoretical counterparts. The infinite number of degrees 
of freedom of the field theoretic systems opens the possibility that the lower 
boundary of the classically allowed region approaches the v = axis only 
asymptotically, or even that it is bounded below by some non-vanishing min- 
imum v (cfr. Ref. |TTJ for a numerical study of classically allowed transitions 



in the Si7(2)-Higgs system). 

3 Quantum-Mechanical Solution 

For the study of the full quantum system it is convenient to use the basis 
formed by the tensor product of the center-of-mass coordinate basis and the 
oscillator excitation number basis: \X) <S> \n). In this basis the state of the 
system is represented by the multi-component wavefunction 

MX) = ((X\®(n\)\y) (11) 

and the time independent Schrodinger equation reads 

- ^QxP + { n + l)^n(X) + £ V n , n ,(X)1> n ,(X) = E^ n {X) (12) 



where 

V n , n ,(X) = {n\\e-^ x+ ^\n') (13) 

In the asymptotic region (large \X\) the interaction terms are negligible, and 
the solution takes the form 

lim VnPO = t±e lk " x + r±e~ lknX (14) 

with 



k n = JE-(n+~)uj (15) 

When n > E — u/2, k n becomes imaginary; we fix the continuation by 
defining 

k n = i^n + ^uj-E (16) 

In order to calculate the transmission probability, we will look for a so- 
lution with 

and 

r+ = (18) 

(This corresponds to an incoming system in its ground state. The general- 
ization to a process with an incoming excited state is straightforward.) The 
inhomogeneous boundary conditions (|17D, (plf) fix the solution completely 
and the transmission probability is then given by 

n<E/u-l/2 K ° 

Numerical methods may be used to calculate the solution satisfying the 
boundary conditions (0), QT7|), and therefore also %. In the rest of this 
section we outline our computational procedure and describe the result. 

To solve the Schrodinger equation numerically we must discretize and 
truncate the system to leave a finite, albeit very large, number of unknowns. 
We accomplish this by replacing the continuum X-axis with a discrete and 
finite set of equally spaced vertices 



Xi = ia 



(20) 



where a denotes the lattice spacing and —N x < % < N x . The truncation in 
oscillator space is perfomed by restricting n < N a . 

To keep our notation concise, we will omit the oscillator indices, using 
implicit vector and matrix notation, and will use subscripts for the locations 
along the X axis. Thus, for instance, the expression J2 n ' Vn,n'('ia)'4>n'{'ici) will 
be simply written as V^. It is also convenient to rewrite the continuum 
equation in the form 

d 2 i)(X) 



dx 2 

where the matrix A(X) is given by 



A(X)^(X) 



(21) 



E 



8n,n' + V n>n \X) 



(22) 



The discretization of Eq. ( plf ) could be accomplished in a straightforward 
manner by using the central difference approximation of the second derivative 
with respect to X 



d 2 ^(X)/dX 2 



ip{X + a) + ip(X - a) - 2ip(X) 



+ 0{a 2 ) 



(23) 



This would give the equations 



Aii/Ji 



(24) 



with an error 0(a 4 ). 

We have actually used the more sophisticated Numerov-Cowling algo- 
rithm, which allows us to improve the accuracy of the discretization by two 
powers of a. From the Taylor series expansion of ip{X) we immediately find 



V'i+i + ^ 



i-1 



d 2 i) 

dX 2 



a 2 + 



ay 
dx* 



— + 0(a b 
.12 v 



(25) 



Using Eq. ©, (<9 4 VV<9X 4 )(a 4 /12) can be written as (d 2 Aip /dX 2 )(a A / '12); 
this can be in turn approximated by (A i+ iip i+1 + Ai-x^i-i ~ 2 A^i) (a 2 / '12) 
with the same level of accuracy. We are thus finally led to the following 
discretization of eq.(BTI): 



a 5(2 (i 

- 2ipi + ^i-i = — A i+1 ip i+1 + —Aiipi + —A i _ 1 'ip i _ 1 (26) 
12 o 12 



which entails an error of order a 6 . 

Before proceeding further, we turn for a moment to the calculation of 
the matrix elements V n ^ n >(X), which are needed for solving the Schrodinger 
equation. These can be calculated very efficiently by means of a recursion 
procedure. V n>n i{X) is given by 

V n , n ,{X) = — 2 KM (27) 
F 

where \v n ) denotes the state 

\ Vn ) = e~^ 2{x+y ^\n) (28) 
It is convenient to use the y coordinate representation, writing 

|^) = e-^^) 2 ^£(^) 1/4 e-^ 2 (29) 

with 



Id [uJ 



2udy V 2 

+ 1 d feu i \ 

a = — 7=3- + \hv 30 

V%jdy V 2 v ' 



The first exponential in Eq. (|29f) can be brought to the right, using 



e -b 2 (x + yf a i = ( flt _ ^£={X + y))e-^ x+ ^ (31) 



This gives 



W = -4=(-) v V - S=(x + y) )\-w«r-w 



The exponential on the r.h.s. may be written as a constant times the ground 
state wavefunction of a shifted oscillator with frequency Q: 

e -h 2 (x + yf-W = Q 1/A e -^Q^ e -^ (33) 

with 

Q 2 

tt = cu + j (34) 



* = V + £x (35) 

Re-expressing everything in terms of the creation and annihilation operators 
for this new oscillator 

1 d [tt 



we find 



, + 1 d n 

b ] = —- + J-Z 36 

s/mdz V 2 K ' 



\v n ) = ^Q lf V^ x2 + (3b + y) n \0) b (37) 



with 



(3 



g 2 



2VujQ 



It is now straightforward to calculate the components of of \v n ) in the h- 
oscillator basis, and therefore also the inner products (v n \v n >), by numerical 
iteration. 

In order to solve the discretized Schrodinger equation, we must calculate 
^ n (with the oscillator index explicit) for — N x < i < N x , < n < N Q . 
This amounts to 2(N X + l)N a complex variables. These satisfy Eq. (|2"6|), for 
N x + 1 < % < N x — 1, for a total of 2(N X — l)N a complex conditions. In 
addition, the boundary Eq. flT7D, (|TB|) give us the 2N a complex conditions^ 

if>- Nmfl = e^V-JV.+i.o + (1 - e tkoa ) 
^- Nx ,n = e tkna ip_ Nx+ltn n>0 

Tp Nx ,n = e^>^-l,n (39) 



2 In these equations we can use either the k n given by the continuum dispersion formula 
Eq. (|l5|) or those given by the dispersion formula that follows from Eq. (pl|). Given the 
high accuracy of the Numerov-Cowling discretization, the two options produce practically 
indistinguishable results. 



Altogether, we thus have a system of 2(N X + 1)N complex, non-homogeneous 
linear equations, which is precisely the number needed to determine all of 
the unknowns. In order to insure good accuracy of the solution, we found 
it prudent to use cut-off values as large as N x = 4096 and N Q = 400. With 
such numbers it would be impossible to tackle the system by brute force 
using a general purpose solver: this corresponds to a system of over 3 million 
complex equations, which could not possibly be solved by a general purpose 
procedure. However, we can take advantage of the special form of Eq. (p6|) 
to implement an efficient solution procedure. Indeed, by inverting a set of 
(N a + 1) x (iV + 1) matrices, which is computationally feasible, Eq. fl2~l~|) can 
be recast in the form 

V>i = Liipi-x + RiA+i (40) 

where Lj and Ri are again (N Q + 1) x (N Q + 1) matrices. The elimination of 
any definite ipi now leads to a system of equations for the remaining variables 
which, with (iV + 1) x (N a + 1) matrix algebra and matrix inversion, can 
be brought to the same form of Eq. (fiOD, with suitably redefined L and 
R matrices. We have used this procedure to progressively eliminate all the 
intermediate variables ipi, i = —N x + 1 . . . N x — 1, ultimately leaving a system 
of equations linearly relating all ipi to ip-N x and ipN x - (Loosely speaking, the 
procedure can be considered the implementation of a Green's function for 
our discretized system of equations.) In particular, ip_ Nx+ i and ipN x -i are 
thus given as linear combinations of ip-N x and ipN x - Substituting these linear 
combinations in Eq. (|39|) we now obtain a system of 2N Q + 2 complex, linear, 
non-homogeneous equations for the 2N Q + 2 complex variables i/j-n x , V^, 
which can be easily solved numerically. As a final remark, we observe that the 
solution procedure outlined above only requires manipulation of real matrices 
for the elimination of the intermediate variables, which entails a substantial 
saving of memory and processor time. Moreover, we can also take advantage 
of the obvious symmetry under reflection of the X-axis to further halve the 
computational costs. 

For our numerical calculations we have used uj = 0.5 for the oscillator 
constant. We have found this value a good middle ground betweeen the 
extremes of very tight and very loose oscillator coupling, where the novel 
features introduced by the internal degree of freedom become less evident. 
Also, apart from some calculations where we varied parameters to study 
the effects of the discretization, we have used a cut-off = 2048 and a 
lattice spacing a = O.OSy^. Insofar as N Q is concerned, we insured that its 
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0.37974509 


0.40829435 


-0.0073 


1.80 


0.34839619 


0.40164508 


0.42390234 


0.43740276 


0.45411756 


0.46898820 


-0.0028 


1.84 


0.52777191 


0.53035499 


0.52913755 


0.52842410 


0.52816695 


0.52887867 


-0.0001 



Table 1: Results for the transmission probability. 



value is large enough that the cut-off energy (N + l/2)u exceeds the barrier 
height by at least a factor of two. We have also checked that the highest 
modes are essentially uncoupled. Specifically, we have used N Q = 400 for 
g 2 = 0.01 and g 2 = 0.02 and N a = 200 for all other values of g 2 (namely 
g 2 = 0.03, 0.04, 0.06 and 0.09). We have solved the Schrodinger equation 
for values of g 2 E ranging between 1 and 2 in steps of 0.02. A slightly thinned 
out compilation of our data is presented in Table 1. 

Our results are also illustrated in Fig. 0. The approach to the classical 
limit (a step function at g 2 E = e (0)) is evident. In Fig. |3| we plot the 
logarithm of the transmission probability as function of 1/g 2 for 8 values 
of g 2 E equally spaced between 1.1 and 1.8. The conjecture following from 
the semiclassical treatment is that for small g 2 the logarithm of the trans- 
mission probability should exhibit the linear behavior in 1/g 2 at fixed g 2 E 
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Figure 3: Logarithm of the transmission probability as function of 1/g 2 . 
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a 


N x 


% 


200 


0.05 


1024 


0.18435 


200 


0.03 


2048 


0.18727 


250 


0.03 


2048 


0.18728 


200 


0.025 


2048 


0.18750 


200 


0.025 


4096 


0.18750 


250 


0.025 


2048 


0.18750 


250 


0.02 


2048 


0.18759 



Table 2: Check of discretization effects. 

(cfr. Eq. (§) 

log To = logC (g 2 E) - \f ( 9 2 E) (41) 

9 

This is well supported by the data in Fig. |3|. We have used the slope of the last 
segment (i.e. the segment corresponding to the two largest 1/g 2 ) for which 
we have significant data to derive the values of F reproduced in Table 1. 
In the next section we will compare them to the results of a semiclassical 
calculation. 

Finally, we performed several checks to estimate the accuracy of our nu- 
merical calculations. For all solutions we verified the degree to which the 
unitarity constraint 

E r(l^l 2 + KI 2 ) = i ( 42 ) 

n<E/w-l/2 

was satisfied. We found this equation fulfilled with an error ranging from 
10 -6 to 10~ 5 . One might object that some of the data for T in Table 1 are 
much smaller than this error. This does not necessarily invalidate them, as 
our use of a discretized Green function may capture the correct exponential 
decays with relative, rather than absolute errors, of the above order of mag- 
nitude. The regularity in even the smallest entries in the Table 1 supports 
this argument. In any event, even discarding all values of T < 10~ 5 one 
would still be left with a rich sample of data verifying Eq. flS] ). 

We have checked the effects of the cut-offs and of the finiteness of the 
lattice spacing by repeating the calculation (for g 2 = 0.03, g 2 E = 1.7) with 
different values of N x , N Q and a. The results are reproduced in Table 2 which 
indicates that, apart from the case of a substantial increase in a, the relative 
errors due to the discretization are of order 10~ 3 . 



An alternative approach to the calculation of % consists in solving the 
time-dependent Schrodinger equation. One can simulate then the collision 
of a wave packet against the barrier and measure directly the trasmission 
probability. It is of course crucial to implement a solution scheme which 
preserves the unitarity of the evolution (up to numerical round-off errors). 
We did follow this approach in some earlier calculations, using a split opera- 
tor technique to achieve a unitary evolution. We obtained results consistent 
with our later calculations based on the time independent Schrodinger equa- 
tion. However solving the time dependent Schrodinger equation proved much 
more (CP) time consuming than solving the time independent one and so 
we abandoned the former method in favor of the technique described in this 
section. 



4 The Semi- Classical Formalism 

We begin this Section with the derivation of the semi-classical procedure for 
calculating the exponent F(g 2 E, g 2 n) of the transmission probability from 
the n-th excited state at total energy E, Eq. (|I0D . Consider an incoming 
state of the form 

\E,n) s = J dP'$p !5 (P')\P',n> (43) 

where 

P = ^2(E - un) , (44) 
\P\ n >= -L / dXe lP ' x \X) ® \n) (45) 

is a simultaneous eigenstate of the center of mass momentum and oscillator 
number, and $p i( j(P') is a momentum space wavefunction with the following 
properties: 

• it is sharply peaked for P 1 « P, with a width of order 5; 

• it corresponds to an X-space wave packet which has support only for 
X <C 0, well outside of the range of the potential. 

With these definitions, the transmission probability is given by 

roc poo 

T n (£) = lim Hm / dX f dy f \(X f ,y f \e- lH ^'^\E } n} s \ 2 (46) 

d^O tf—U— >oc Jo J-oo 



This motivates us to calculate the matrix element 



A(X f ,y f ,P,n) = (X f ,y f \e-* H <tt-V\P,n) (47) 

Position-eigenstate matrix elements may be evaluated in terms of a path 
integral involving the classical action: 

(XfMle^t-VlXiM) = C J[dX}[dy]e lS (48) 

where C is a normalization constant, and the integration is over paths sat- 
isfying X(ti) = Xi, y(ti) = yi and X(t f ) = X f , y(t f ) = y f . The amplitude 
P7| ) is the convolution of the path integral ( |4"8"[ ) with the eigenfunctions of 
the center-of-mass momentum and oscillator excitation number, e tPX and 
(y\n), respectively. (y\n) is conveniently represented in terms of an integral 
over coherent state variables z and z. In this way we obtain 

A(X f , y f , P, n) = -±= / dX % d Vi e* px * f ^ e^^= 

e -$*-M-^*vi(X f ,y f \e- %H <t'- ti ')\Xi,y i ) (49) 

The main idea, adapted from the method of Ref. [0, |[, is to set E = e/g 2 , 
n = v 1 1 g 2 and take the limit g —>■ while holding e, v fixed. Indeed, by 
rescaling the integration variables, we are then able to recast the matrix 
element in the form: 

A = cjdX l dy l J^ J [dX] [dy] e'^ T (50) 

with 

1 i 1 

T = — iS + -z 2 + T^Vi — \ 2tuzyi + zz — ipXi — v In z + 2^( m v ~~ ^) 

where p = gP, S is the classical action, and Stirling's approximation has 
been used for the factorial^. 

This form for the matrix element A is now suitable for a semi classical 
analysis at small g: we find stationary points of T, and evaluate the integral 

3 In order to keep our notation simple, we have used the same symbols (i.e. Xi, yi) for 
all rescaled integration variables. Note that Xf and y/ must also be eventually integrated 
upon, cfr. Eq. (f46|), and are rescaled as well. 



in a gaussian approximation about such points. The stationarity conditions 
(obtained by varying X(t),y(t),Xi, z, z, and yi) are: 
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SX(t) 5y(t) 

dX 
~dt 



0, t^u,t f 
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(52) 
(53) 
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2uj z 
2uj z 



(54) 



As expected, Eq. Q52|) is the classical equation of motion for this sys- 
tem, while Eq. (|53|) and (|54]) imply that the initial classical state has the 
(rescaled) center-of-mass momentum p and oscillator excitation number u, 
so that e = p 2 /2 + up. Therefore in the classically forbidden region of the 
e — v plane there will be no real solution where the system goes over the 
barrier. Nevertheless there may be complex solutions. We expect that the 
integral is dominated by stationary points, even if these points lie outside 
the domain of integration. Hence we will seek solutions that may involve 
complex values for the integration variables^. In searching for such solutions 
we must remember that we are performing an analytic continuation of the 
integration variables; in general we will run into singularities in the com- 
plex t-plane. To deal with this problem we note that the time contour, the 
real time axis, can be distorted into the complex plane without changing the 
path integral, provided we keep the time contour end points (ti,tf) fixedF]. 
Thus our strategy will be to search for complex solutions to Eq. ([54]) along 
a complex time contour ABCDE as shown in Fig. ^. 

The matrix element ( |4"T| ) in this approximation becomes 



.4 



(55) 



4 In general this allows values of z and z such that z ^ z* . 
5 The evolution operator may be written exp[— iH (tf — ti)] 



(tf — ti). This argument applies even for complex dtj 



Y[j exp[— iHdtj] provided 



Imt 




Figure 4: Contours in the complex time plane used to find the saddle point 
solutions. 



where the correction c determines a pre-exponential factor in the semiclassical 
limit, i.e. lim g _ t0 g 2 c = 0. The function T is equal to the right hand side of 
Eq. fl5"T|) evaluated at the solution to Eq. ( |52| - p3D . Using these equations to 
eliminate z we write 



1 1 
-iSq — — v In z + —v In v 



(56) 



where 



dt 



2 / dt 



t=t f 



1 dy 



t=t f 



(57) 

Note that the quantity So is insensitive to the value of ti, provided the solution 
to Eq. ( |52"D is in the asymptotic region near time ti. 

Equations (|53D, (0) and (^) involve the quantities defined at large neg- 
ative time on the real time axis (region A in Fig. ^). It is convenient to 
formulate the boundary conditions at large negative Re t on the part BC of 
the contour, where 



t = t' + -T , t' = real 
2 ' 



-oo 



Since for the moment we consider the asymptotic past, we may ignore the 
potential V and write the solution at large negative t' as follows, 



y(t) 



(58) 



'2uj 

X(t)=X Q +pt' (59) 

The three parameters of the solution, Xq, u and v, which are in general 
complex, are related to the quantities entering Eq. (|53D , (0), (|56]) in an 
obvious way, 



ue 2 



uiT 



ze 



ve 2 



uT 



Xi=X - -pT + pti 

The condition 

ZZ = UV = V 

tells us that the phase of u is opposite to that of v, and we may parametrize 
them as 

v = eV (61) 



(60) 



So far we have not specified a value for the parameter T; we have argued that 
our result is independent of this parameter, provided we avoid singularities 
in the complex plane. Since variation of T changes the value of Im Xq, we 
can adjust T such that X is real in the region B: Im Xq = 0. 

The transition probability is given by the absolute value of the matrix 
element squared; that is, in terms of twice the real part of T. Using the 
above relations between u and v we can write the transmission probability 
as exp (—F/g 2 ) where 

F = 21m S - eT - v9 (62) 

The resulting values of T and 9 depend on v and e. However, we may treat 
T and 6 as independent parameters instead, so that the boundary conditions 
are formulated in a simple way in the asymptotic past on the part BC of the 
contour: 

(i) X(t') and X(t') are real at B. 

(ii) positive and negative frequency parts of the oscillator solution fl58 ) 
are related by Eq. (|6"lD at B. 

At given T and 9, the initial center-of-mass momentum and excitation 
number (and hence the total energy) are to be found from Eq. Q53| ) and 
It is straightforward to check that 

d(2lmS (T,9)) = 
dT 6 

d (21m S (T,9)) 



v 



39 

so that T and 9 are Legendre conjugate to e and v. It is worth noting also 
that Eq. ([Jl]) at 9 ^ in fact requires the solution to be complex in the 
region B of the contour. 

We are interested in the total probability for transmission; thus we should 
integrate our probability over all values of jjf and over positive values of Xf. 
This final integral may also be done using the saddle point approximation. 
The saddle point condition is simply that 

(iii) the solution X(t) and y(t) should be real along the D —> E part of 
the contour. 

At given T and 9, the classical equations of motion and the boundary 
conditions (i), (ii) and (iii) are sufficient to specify the complex solution up 
to time translations along the real axis. Finding the solutions is still a non- 
trivial computational task. To simplify this task we start from a sub-class 



of solutions with 8 = 0, whose numerical determination is easier, and then 
deform these solutions to 9 ^ 0. 

For the solutions with 9 = the X and y coordinates are analytic and 
real along the entire contour BCDE of Fig. |j. From the Cauchy-Riemann 
conditions it follows that the velocities X and y vanish at C and D. The 
motion along the imaginary time axis can be reformulated in terms of r = 
Im t and a "Euclidean" Lagrangian 



(63) 
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Figure 5: Motion in the "periodic instanton" solutions. 



The equations of motion 

-(X + y)e-^ x ^ 2 

u 2 y - {X + y)e-^ x+y)2 (64) 

describe the evolution of the system along the imaginary time axis. We look 
for periodic solutions where dX/dr and dy/dr vanish at the turning points 
t = and t = T/2, T being the period of the motion. These solutions are 
analogous to the Euclidean solutions that are commonly used to describe 



d 2 X 
~d^ 

dr 2 



motion through a barrier in the semiclassical treatment of tunneling with 
a single degree of freedom. In this latter case finding periodic solutions is 
straightforward: one need only integrate the equations of motion with in- 
verted potential. The situation with several degrees of freedom is not so 
simple. Indeed, the continuation to imaginary time not only inverts the po- 
tential barrier, which now becomes a potential well, but also changes the 
harmonic restoring force into a linearly increasing repulsive force. This force 
makes the system unstable, requiring careful adjustment of the values of X 
and y at the turning points to obtain a periodic solution. The resulting mo- 
tion is similar to the "periodic instanton" solutions that appear in topology 
changing transitions in quantum field theory [12]. There too, all of the field 
oscillator degrees of freedom become repulsive in the Euclidean motion and 
the solutions are unstable: a small perturbation of the field profile at one 
of the turning points grows exponentially in the subsequent evolution. The 
motion in the "periodic instanton" solutions of our model is illustrated in 
Fig. |5|. At the turning points (Pi, P5) particle 1 is attracted towards the bot- 
tom of the potential well but repelled by particle 2, which is located between 
particle 1 and the bottom of the potential. Both particles accelerate towards 
the bottom (particle 2 because of the repulsive force exerted by particle 1). 
The balance of forces, however, is such that particle 1 moves faster than 
particle 2, reducing the interparticle distance (P 2 , P4) and, correspondingly, 
the repulsive force, until it overtakes particle 2 precisely when both particles 
transit through the bottom of the potential (P3). 

Finding the periodic instanton solutions of our model is rather easy. We 
divide the interval < r < T/4 into N subintervals of width At = T /(AN) 
and we denote by Xi, yi the values taken by X and y at r = iAr. We 
discretize the Euclidean action Se = j ' Le<It and look for a minimum of Se 
with respect to the variables X^y^ i = . . . iV — 1, while X N and y^ are 
kept fixed at zero. Since the Euclidean action is bounded from below, the 
algorithm of conjugate gradients converges rapidly to the correct solution. 
The values Xo,yo can then be used as initial data for the integration of the 
equations of motion along the real time axis (from D to E in Fig. |J) until both 
particles are far out of range of the potential. In this manner one can find 
the asymptotic oscillator number of the solutions (in the periodic instanton 
solutions initial and final oscillator numbers are of course identical). The 
periodic instanton solutions span the one-dimensional subspace denoted by 
the thick line in the v-e plot of Fig. 0. 
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Figure 6: The curve spanned by the periodic instanton solutions. 



Starting from the periodic instanton solutions we find solutions with 6 > 
and reduced incoming oscillator number v by a deformation procedure. In 
order to solve the equations of motion numerically we subdivide the contour 
BCD of Fig. |] into iV subintervals separated by vertices labeled by an index 
i = . . . N. We denote by 5ti the time interval from vertex i to vertex 
i + 1. 5ti will be real along the BC portion of the contour, after which 5ti 
should be negative imaginary. It is convenient, however, to place the very 
last subinterval along the real time axis. Thus we place the vertex iV — 1 at 
the origin, and the vertex N at the point t = 5tN-i = real. The discretized 
action is 
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N-l 

E 

i=0 
— UJ 



{X 
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x 4 
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2SU 

\(X l+1 +y l+1 ) 2 _|_ e -\(X x +y,f 



(65) 



This expression is quite general and valid for any contour of integration in 
the complex time plane. 



It is convenient to use Zij, j 



,4, to denote the four variables 



Re Xi, Im Xi, Re yi, Im j/j. The equations of motion are given by 

dS 



dz. 







1...N-1 



'■j 



These amount to 4N — 4 conditions for the 4N + 4 unknowns z, 
solution must also satisfy the boundary conditions 



h3 ■ 



2/o + 2/i — 
[cfr. Eq. (H), (0)] and 
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Im X 



N-l 



Im X/v = Im yN-i = Im t/iv = . 



(66) 
The 

(67) 
(68) 



In addition, we remove the invariance under time translation by demanding 
that X takes a fixed real value 



^0 



(69) 



The precise value of c is not relevant. The only important criterion that 
c must satisfy is that the imaginary time axis falls between the expected 



singular points of the solution. The value of c can be readjusted, if necessary, 
so that the point in the complex time plane where ReX = belongs to the 
CD part of the contour. 

Equations (|6?D-(j69|) provide the required 8 additional conditions on the 



(70) 



variables 2j,-. We will write these equations as 







Starting from the periodic instanton solutions and evolving them further 
from C to B we obtain an initial class of solutions to Eq. (|66|), (|T0]) with 
= 0. If we perform a small change of either 6 or T, the two parameters 
which indirectly determine v and e, the field configuration Zij will no longer 
satisfy the equations of motion. We seek a correction Szij such that Zij + 5z it j 
obey the equations of motion with new values for 9 and/or T: 



OS 



dzij 







(71) 



\-Sz 



B k ( Zijj + Szij) = (72) 

If the deformation of the original solution is not too large, Eq. (|7T|), (|72|) can 
be solved by the Newton-Raphson method: we expand to first order in 5z 
and solve the linearized equations 
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j, dzi>j> 



QZ%i j> 



(73) 



(74) 



This procedure is repeated until it converges to a solution. 

In our calculations we typically took N = 2048 and used the following 
computational strategy. Equation (|73|) with a definite index i only couples 
the variables SZi'j with i! = i — 1, i, i + It is then possible to use an elimi- 
nation procedure similar to the one outlined in Sect. [3] (see Eq. (^0|) and con- 
siderations that follow) and express all variables Szij in terms of Szqj, 5zn,j- 
(In practice this can be done maintaining complex variables notation, which 
simplifies the arithmetic. One must work with the explicit real and imaginary 
parts of the variables only at the next stage of the calculation.) Finally Szoj, 
Sznj, and 8z 1 j, 5zjv-ij which, by virtue of the elimination procedure are 



now expressed as linear functions of Sz j, Sz^j, are inserted into Eq. ([73]). 
These equations thus become a system of 8 real, linear, non-homogeneous 
equations in the 8 real variables Szqj, Sznj, that can be straightforwardly 
solved. We are then able to start from 8 = (periodic instanton solution) 
and gradually increase 9 to a very large value, which makes the incoming 
oscillator number v effectively zero. At the same time we gradually reduce 
the value of T, which has the effect of increasing the energy e. It is important 
to check that the solutions correspond indeed to tunneling processes, namely 
that in the further evolution along the positive real time axis the center of 
mass coordinate X goes to +oo. We found this to be the case up to e ~ 1.1. 
At that point, though, the Newton-Raphson method develops an instability 
and, when convergence is eventually reached, further evolution along the real 
time axis shows a bounce from the barrier with X — > — oo. We attribute this 
difficulty to the proximity of solutions with X — > +oo and X — > — oo, with 
possible bifurcation points. In order to avoid falling into a solution without 
tunneling, we continue a tunneling solution to a positive real value of t along 
a contour extending into Im t < 0, as illustrated by the dashed line in Fig. |j. 
At the final value of t (point E in the graph of Fig. ^) the system is far into 
the positive X domain and we can further increase e without running into 
any instability. 

We illustrate in Fig. [7] a typical tunneling solution in the complex time 
plane. The figure displays the center of mass coordinate X as function of 
Re t, Im t (we inverted the direction of the Im t axis for a better perspective). 
The height of the surface gives the value of Re X, while the phase of the com- 
plex variable X is coded by color (red for real negative, blue for real positive, 
with the other values of the complex phase arranged in rainbow pattern — the 
color will appear as gray-scale in a black and white printout). The contour of 
integration of the equations of motion is indicated by a line of different color 
drawn on the surface. The continuation of X to the entire complex plane 
has been obtained by starting form the values along the imaginary time axis 
and integrating the equations of motion outward with the leapfrog algorithm. 
The singularities in the solution are quite apparent from Fig. |7J. We found it 
noteworthy that one can determine the singularity structure of the solutions 
numerically, since ordinarily one would expect numerical integration methods 
to fail in the presence of a singularity. The integration algorithm becomes 
unstable and diverges as one approaches a singularity. However it is possible 
to exhibit the singularity structure by numerically integrating the solution 
along closed contours around the singularities. Integration of the equations 



Figure 7: Tunneling solution to the equations of motion: center-of-mass 
coordinate as function of complex time. The singularities (branch points) 
are labeled by crosses. 



of motion by the leapfrog algorithm along a closed contour entirely contained 
within a domain of analyticity produces final values for X and y identical to 
the initial values to high degree of numerical accuracy, equal to the expected 
discretization error 0((5t) 2 ) of the algorithm, whereas an enclosed singularity 
is clearly present when the initial and final values of X and y are different. 

We reproduce in Table ^| the results of our semiclassical calculation. The 
data correspond to 9 = 13. In Figure § we present a comparison of the 
results for the exponent F in the transmission probability (cfr. Eq. ([!])) ob- 
tained with the full quantum-mechanical calculation (x) and with the semi- 
classical technique (solid line). In the quantum-mechanical calculation we 
extracted F from the slope of the last segment in the graph of log % versus 
1/g 2 at given energy for which we had meaningful data (see Fig. |3|). As a 
consequence, the line defined by the crosses in Fig. |8] exhibits some small 
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Table 3: Results of the semiclassical analysis. 
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Figure 8: Comparison of the quantum mechanical and semiclassical results. 



discontinuities. We take the magnitude of these discontinuities as an indi- 
cation of the systematic errors in the quantum-mechanical calculation due 
to the neglect of perturbative 0(g 2 ) and higher order effects. Within these 
errors, the agreement between the results of the full quantum-mechanical 
calculation and of the semiclassical calculation is excellent. 

5 Conclusions 

Our results validate, in the context of a model calculation, the scaling for- 
mula of Eq. ([J), @, the applicability of the method of Ref. |], || and the 
assumption that the ground state transition probability can be obtained as 
the limit of a more general transition probability from a coherent initial state. 

At the same time our investigation has brought to light interesting prop- 
erties of the analytic continuation of classical solutions to complex time and 
complex phase space. While the extension of classical motion to the com- 
plex time domain has long formed the mainstay of semiclassical calculations 
of tunneling, we believe that our specific application shows novel features 
of the analytically continued solutions intimately connected to the presence 
of several degrees of freedom. Of particular relevance we find that one can 
obtain information on the singularity structure of the solutions by numerical 
techniques. 

With the qualification that a field has an infinite number of degrees of 
freedom while our model has only two, our results bode well for the applica- 
tion of the technique of Ref. |], |f to field theoretical processes. Hopefully, 
they will also open the path to new, imaginative applications of semiclassical 
methods in other challenging quantum-mechanical problems. 
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